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Abstract 

We analyze two- and three-link planar snake-like locomotion and optimize the motion for efficiency. 
The locomoting system consists of two or three identical inextensible links connected via hinge joints, and 
the angles between the links are actuated as prescribed periodic functions of time. An essential feature 
of snake locomotion is frictional anisotropy: the forward, backward and transverse coefficients of friction 
are different. The dynamics are studied analytically and numerically for small and large amplitudes of 
the internal angles. Efficiency is defined as the ratio between distance traveled and the energy expended 
within one period, i.e. the inverse of the cost of locomotion. The optimal set of coefficients of friction to 
maximize efficiency consists of a large backward coefficient of friction and a small transverse coefficient 
of friction, compared to the forward coefficient of friction. For the two-link case with a symmetrical 
motion, efficiency is maximized when the internal angle amplitude is approximately 7r/2, for transverse 
coefficient sufficiently large. For the three-link case, the efficiency-maximizing paths are triangles in the 
parameter space of internal angles. 

1 Introduction 

Snake locomotion has long been a topic that fascinates researchers, and has recently received a renewed 
wave of interest in the fields of robotics and control, as well as in organismal biology. Snakes are familiar 
organisms, but as limbless animals, their locomotion has special features [7j- Terrestrial snakes move using 
friction between the ground and their belly scales, which have anisotropic frictional properties [5^. It has 
been proposed that the cost of transport (energetic efficiency) for snake slithering is no greater than that of 
limbed animals [251 [?7]. 

Some works on modeling snake locomotion are oriented towards wheeled snake robots [31 21 |S1 IHl 
[TTl [T^ . These models are typically concerned with motion planning, and assume that the transverse 
coefficient of friction is high enough to prevent transverse motion, while the forward and backward coefficients 
of friction are the same. These models work well for wheeled robots and provide valuable insights into the 
locomotion of biological snakes. In experiments, Hu et al. [12^ measured the frictional anisotropy of juvenile 
milk snakes, and found that the forward, backward and transverse coefficients are different but similar in 
magnitude. They also studied the effects of the snake's active modulation of its weight distribution on 
the ground. Hu & Shelley [13^ analyzed the so-called "lateral undulation" motion modeled as a family of 
sinusoidal traveling-wave shapes and calculated the dependence of speed and efficiency on amplitude and 
wave length of the kinematics as well as coefficients of friction. 

In this work, we adopt the model of ^12] , but for bodies composed of two and three rigid links. Linked 
bodies are fundamental for robotic sliding systems. By specializing to bodies with two and three links, 
we consider the simplest such systems, which nonetheless have nontrivial behaviors. Two- and three-link 
locomoting bodies have been considered previously as swimmers at zero Reynolds number (Stokes flow). 
Purcell [22) described the physics of swimming in Stokes flow, and stated the Scallop Theorem: in Stokes 
flow, net locomotion is not possible if a swimmer deforms in a way that is invariant under the reversal of 
time [5 . Such is the case for periodic motions of a two- link body ("scallop"). He then proposed a three-link 
swimmer that moves only one link at a time, in a non-reciprocal motion that results in net locomotion. 
Subsequent studies have calculated efficiency-optimizing motions for Purcell's swimmer and similar systems. 
Becker et al. j2j calculated efficiency-optimizing stroke amplitudes for Purcell's swimmer, and considered 
different length ratios of the three links. Tam & Hosoi extended the optimization to arbitrary kinematics 



(allowing both internal angles to change simultaneously) and arbitrary slenderness ratios. They found an 
optimal path in the parameter space of internal angles using a Fourier series representation, and showed 
that the high frequency modes are subdominant to the low frequency modes. Avron & Raz [1] developed 
a qualitative geometric approach by focusing on the curvature of the local connection matrix to study the 
same system. Hatton & Choset [TUl further developed this technique, and suggested a systematic way of 
choosing the best body-fixed frame to approximate the inertial frame displacement while accounting for the 
overall rotation. They calculated optimal motions for other systems such as a three-link fish swimming in 
infinite Reynolds number (potential flow) , which admits a similar formulation. Kanso et al. [15] and Melli et 
al. |18) gave a geometric formulation for swimming in a potential flow and calculated optimal strokes. Jing 
& Kanso [13] used this formulation to study the effects of elasticity and body configuration on the stability 
of passive locomotion. 

Here we study the locomotion of two- and three-link bodies not in a fluid but instead on a planar 
surface with sliding friction. Unlike the fluid studies, the dynamical equations are nonlinear with respect to 
body velocities, so the Scallop theorem and many other results no longer apply. Like the aforementioned 
studies, we focus on finding motions which optimize the efficiency of locomotion. The structure of the paper 
is as follows. In Section [2] we formulate the model for a two-link "snake," nondimensionalize the system 
and derive the equations of motion. In Section [3] we use analysis and computation to optimize the two-link 
model with respect to kinematic parameters and coefficients of friction, for both small- and large-amplitude 
actuations. In Section |4] we analyze the three-link snake model and compute optimal kinematics of the 
relative angles at one realistic set of coefficients of friction using a Fourier series representation. In Section [5] 
we summarize our work and suggest directions for future study. 

2 Problem setup for two-link model 




Figure 1: Two-link snake model; see text for description. 

The snake is modeled in 2D as two identical inextensible line segments (links) connected via a hinge 
joint as depicted in Figurejl] The total length of the snake is L, and for each link is L/2. The snake shape is 
parametrized by the angle 9r between the tail link (left) and the head link (right), with the positive direction 
of rotation being counterclockwise. Denoting the snake's mass per unit length as p, the total mass is m = pL. 
We denote s as the arc length between any point of the snake and the tip of its tail, so < s < L. The tail 
tip, hinge joint and head of the snake correspond to s = 0,i/2 and L respectively. 

Motion of the two-link snake can be observed both in an inertial frame {Gj., e.y} with origin at a fixed 
point O, or in a body-fixed frame {b^., bj^} with origin at the center of mass C. The unit vector b^, is parallel 
to the line connecting the links' centers, and hy is rotated by 90 degrees. The position of C in the inertial 
frame is Xc = (xc, j/c), and the orientation 9c is the angle from to b^.. In the inertial frame, the position of 
an arbitrary material point on the snake is denoted as x = {x,y), and 6 is the angle between the tangent to 
the snake at a given point and e^;. The positive direction of rotation is counterclockwise. In the body- fixed 
frame, the position of the same material point is denoted X — {X,Y), and the tangent angle is Q. For a 
given material point, we define the configuration variable in the inertial frame g = [x , y, 9]'^, and in the 
body-fixed frame G = [X , Y , 6]^. Specifically, for C, we have gc = [xc, Vc , 9cY' and Gc = [0,0, 0]-^. 
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Figure 2: Schematic of the body-frame velocities for the two-link model. Center-of-mass velocity is = [Uc i Vc i f^c]"^, the 
velocity due to rotation about C is $,rot = {—VtcY , ficX , 0]^, where {X , Y) is the position of a material point in body frame, 
and the velocity due to shape change ^shape consists of a horizontal motion of the link centers and rotations about the link 
centers with angular velocities ±8r/2 for the two links (in this example, 8r > and 6r < 0). 



The relation between the configuration variable in both frames is 

/ cos 6'c - 

g = gc + Re^G, Re^=\sinec cos^e 0|, (1) 




where Rg^ is the transformation matrix. For the tail link, i.e. < s < L/2, the configuration in {b2:,bj,} is 
given by 

fxA ( [s - 1/2) cos(^,/2) \ 
Gt- = -(s-l/4)sin(0,/2) , (2) 

\Qtj \ -Or/2 ) 

where the subscript t indicates tail. For the head link, L/2 < s < L, the configuration is 

^(s- l/2)cos(6'^/2)^ 
Gh =\Yh\ = \{s- 3/4) sin(0,/2) I , (3) 





er/2 

where h represents head. 

In the inertial frame, the linear and angular velocities of any point on the snake are given by the 
time derivative of its configuration, that is g. In particular, the velocity of C in the inertial frame is given 
by gc- The velocity of C with respect to the inertial frame can also be expressed in the body frame as 
= \Uc , Vc , ^c]'^ , where ^Ic = 9c- Similarly, the velocity of any material point with respect to the inertial 
frame can be expressed in the body frame as 

^ = 1^1 = shape +Crot, ^shape = ' \ Irot = ( ] , (4) 

where $,shape is due to shape changes and ^rot is due to rotation of the snake about C (see Figure [2|. The 
relation between the velocity in both frames is given by 

g = ReA^ in particular gc = ReAc- (5) 

Perpendicular to the plane of motion, the forces on the snake are gravity and the supporting force 
from the ground, which balance each other. Within the plane of motion, the forces on the snake are external 
friction from the ground and internal forces. Since the coefficients of friction are anisotropic, the frictional 
force is decomposed into components in different directions. In the body frame, we denote the linear velocity 
of an arbitrary material point as = [U , V]'^ , the unit vector tangent to the snake as s = (cos 0, sin 0), 
and the unit normal vector as n = (— sinO, cos 8). Our model for the snake mechanics is essentially the 
same as that in |12j . The Coulomb frictional force density at a given point on the snake is 

f{s,t) ^ pg ^-fitiihn ■ ii)n- fJ.fH{iiin ■ s) + Hb (l - H{iun ■ s)j (^Hn-s)s|, (6) 

where fJ.f, fJ.b and fit are the forward, backward and transverse coefficients of friction respectively, = 
iiin/Miin\\, and H{-) is the Heaviside function, used to distinguish forward and backward friction [12]. Note 
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that the frictional force density depends on the direction of velocity but not the magnitude. In addition to 
the external force, there are also forces internal to the snake, and the internal force density is denoted as 
fin{s,t). The torque density with respect to C due to friction is given by T{s,t) — X-*^ • f, while that due 
to internal force is Tin{s,t) — X-'- • fi„. The internal force and torque densities are due to a system of equal 
and opposite tension and shearing forces acting on adjacent sections of the snake across their interfaces [M) . 
which makes the snake inextensible and enforces its shape. The integrals of internal forces and torques are 
zero 

/•L pL 

/ f™ds = 0, / n„ds = 0. (7) 
Jo Jo 

Now we nondimcnsionalize the variables. We consider actuations 6r (and resulting snake motions) 
that are periodic in time with period T. Variables are nondimensionalizcd by scaling by the total length L, 
period T and mass m = pL. Three important dimensionless numbers for the snake dynamics are 

Here Fr is the Froude number, which can be written as a ratio of snake inertia to the forward frictional force 
acting on it. The other two parameters are friction coefficient ratios. We assume the coefficients of friction 
are uniform along the snake, and define the forward direction as that with the smaller of the tangential 
friction coefficients (if p,f ^ Hb), so /ib > 1 by definition. For real snakes, Fr <C 1, which means that the 
snake's inertia is negligible compared to frictional forces [T21 113j . 

For simplicity, we now drop the tildes with the understanding that all variables are dimensionless in 
the remainder of this work. The dimensionless point-wise frictional force density is given in the body frame 

by 

f = -fJ-tiiun ■ n)n - ^Hiiiin • s) + 1 - i?(|"/„i • s) | {iun ■ s)s. (9) 

It can be transformed into the inertial frame via the transformation matrix Rg^. The governing equations 
are given by the linear and angular balance laws, 

Re^ ^ (f + f„, ) ds = X, , ^ X^ • (f + f„) ds = ^ ^ {Fr x iu„) ds . (10) 

Since the inertia term is negligible, the Froude number is assumed to be zero: Fr — 0. Therefore, substitut- 
ing ([t]) into (10), the governing equations are reduced to the integrals of frictional force and torque densities 



equaling zero, 

»i .1 
fds = 0, / X-L-fds^O. (11) 
'0 Jo 

The above equations give three independent scalar equations for the three unknowns = [Uc , K , ^c]'^ , 
and the solution depends on the parameters fib and fit, as well as the prescribed relative angle 6r{t). Notice 
that (11) are intrinsically nonlinear, and the nonlinearity primarily lies in the form of friction given in (|9|. 



An explicit derivation of the equations of motion is given in Appendix [Xj 

Since f only depends on the direction of the velocity, and no inertia term is present in the equations, 
the only time scale in the problem is T. For a given 9r{t), if T is doubled, the speed of the motion is 
reduced by half, but the snake will trace exactly the same trajectory, as does Purcell's three-link swimmer in 
Stokes flow. This is referred to as the rate independence or time invariance of inertialess systems: if a body 
undergoes a deformation, the trajectory traveled by the body between two different shape configurations 
does not depend on the rate of deformation, but only on the sequence of deformation [2j [5l [22] . On the other 
hand, the Scallop theorem indicates that inverting the shape change sequence corresponds to inverting time 
for the original shape change sequence (note this is not equivalent to time invariance) . A corollary [1] [TU] 
is: in the body-fixed frame, the velocity of the center of mass is proportional to the velocity of the shape 
change. This is known as the kinematic reconstruction equation in the geometric mechanics literature 

^, = A{er)er, (12) 

where 0^ is a shape change vector for systems with multiple degrees of freedom, A{6r) is the local connection 
matrix, which is an n x m matrix that relates an to x 1 shape change velocity 0^ to an n x 1 velocity of C 
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in the body frame ^c- For the two-hnk model, 9r is a scalar and is a 3 x 1 vector, hence A{6r) is a 3 x 1 
vector. Note that A{9r) does not depend on 9r. 

However, equation ( 12 1 does not apply to our snake model. This is due to the anisotropy of coefficients 
of friction, specifically, the Heaviside function in the force equation ^ , that causes the irreversibility of shape 
change. Instead, for our system, the local connection matrix also depends on the direction (or sign) of Or, 
denoted by Sr — sgn(0r), but not the magnitude ||^r ||- If there were no Heaviside function in ([o]), and force 
were only decomposed into tangential and transverse directions, then the system would become very similar to 
the multi-link fish problem in a potential flow, for which it is known that with the added inertia decomposed 



into tangential and transverse directions (12 1 holds jl5j . Consequently, the techniques of analyzing the local 
connection matrix developed in pfl [TO] cannot be directly implemented here since it is based on (12). In 
general, the modified kinematic reconstruction equation for the snake model can be written as 



(13) 



For the two-link model where 9r is a scalar, the above equation is reduced to — A(6'r, Sr)dr- Note when 



relative angle 6r{t) is prescribed, (13 1 is a nonlinear algebraic equation rather than a differential equation. 



and the nonlinearity lies in the form of A(6',., Sr)- The solution of (111 or equivalently ( 13 ) is the velocity of 
C expressed in the body frame, ^c- Without loss of generality, assume the snake starts from the origin with 
zero initial orientation angle, i.e. gc(0) = 0. The configuration in the inertial frame is then 



Jo Jo 



(14) 



which is an iterative integral in time since 9c is only known from the integral 9c = /q ^c^) dt. The distance 
traveled by C during one period observed in the inertial frame is given by 



d = ||x,(r)-x,(o)|| 



(15) 



The work generated by the snake during one period is equal to the energy dissipation due to friction since 
the system is inertialess, i.e. 



W 



1 ■ iun ds dt. 



(16) 



here it is more convenient to express f in the body frame. The efficiency of locomotion is defined as the ratio 
between distance and work, 

"w- ^ (-) 

This efficiency e, after nondimensionalization, is equivalent to the inverse of the cost of locomotion commonly 
seen in the animal locomotion literature [5 . Intuitively speaking, e is analogous to the concept of "milcs- 
per-gallon" . 




48„ 



(a) 6r vs t 



0.5 
(b) Br vs t 



Figure 3: (a) Relative angle 6r prescribed as a triangular wave with period T = 1 and amplitude 9max- (b) Angular velocity 
for relative angle 0r. 
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3 Two-link model analysis 



We first look at an example of two-link snake locomotion. We prescribe the relative angle 6^ as a triangular 
wave with period T = 1 and amplitude Omax, and with 9r{Q) = dmax, i-e. 



0rit) 



,,,(l-4t), 0<t<0.5, 
,,,(-3 + 4i), 0.5<t<l, 



Q<t< 0.5, 



^Oraax, 0.5 <t< 1. 



(18) 



1.3 and fit — 1-7, which are 



Figure [s] depicts 9r and Br within one period. For concreteness, let fib 
taken from the experimental measurement in 13J. The equations of motion (111 are solved numerically 
using the subroutine f solve in MATLAB, which implements a Trust-Region-Dogleg algorithm. When 
Omax = 7i'/2, the trajectory of the center of mass C is shown in Figure Qa) with five snapshots of the snake 
at i = 0, 0.25, 0.5, 0.75 and 1 overlaid (the head of the snake is represented by O), and the orientation 9c as 
a function of time is depicted in Figure |4](b). One can see from the trajectory plot that the displacement 
in the x direction is larger than that in the y direction. And even for the large amplitude of actuation 
0max = T^/"^, the rotation is very small: \\9c\\ < 5 x 10"^ = 0.3 degrees for all time. The distance traveled 
during the period is d = 0.0535, the work is ly = 0.7868, and the efficiency is e = 0.0680. Figure [5] shows 
the velocity of the center of mass gc in the inertial frame. The horizontal velocity Xc is always non-negative. 
The vertical velocity ijc is non-negative during the first half-period, and non-positive during the second half. 
The two half periods nearly cancel out and result in nearly zero net displacement in the y direction as shown 
in Figure |4]( a). The angular velocity 6c alternates between positive and negative during the period, and is 
symmetric about t — 0.5. All components of velocity become when t = 1/4 and t — 3/4, which corresponds 
to 9r — 0. From Figures [4] and [H] one can see the resulting velocities depend nonlincarly on 9r. They can 
only be solved numerically. The nonlinearity arises from both the form of the friction in Q, and the large 
amplitude of actuation, in this case 9max — tt/2. In order to focus on understanding the former nonlinearity, 
we now analyze an actuation with small amplitude. 
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(a) trajectory of C and snapshots 



t 



o.r, 1 

(b) orientation 8c vs. t 



Figure 4: Motion of two-link snake when dr is prescribed as a triangular wave given in | |18| l and dmax = tt/2: (a) Trajectory 
of center of mass C in the inertial frame with snapshots of the snake at t = 0,0.25,0.5,0.75 and 1; (b) Orientation 0c as a 
function of time. The coefRcients of friction are fib = 1-3 and /it = 1.7. 
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Small amplitude analysis For a general kinematics 9r{t), assume sup^ \\0r\\ = e ^ 1, and consequently 
Or, Or ~ 0{e). One can show that Uc, Vc, flc ^ 0{e) as follows. Specifically when the motion of the snake 
is as depicted in Figure 
equations of motion (11| results in 



ire 2j that is 0^ > and Or < 0, we show in Appendix |A| that the integral form of the 



1 

16^ 



fJ-b [ Uc + —OrOr I In 



iUc + OrOr/4: 



Uc 



16 



Or In 



-Or 



AUc - OrOrl^ 



16 



-at {2Uc0r - ^c) 



( 2K + ^jrOl ) - A^b ( UcOr + Y^OrOl ) hi 



4^t(l + 



4Uc 



(1 - /i6)6'r 



9./4 



(19) 



1. 
16 



Or In 



The solutions for equations (|T9|) are given by 
'Uc 



1 



--(l + /3)0. 



Vc 



^c 



1 

32 



-1 



/31n 



(20) 



-^(l + /3)-^/3[/3(l + A'b) + 2M6] In 



where /3 is the solution of the following transcendental equation 



[ipb + l)P{t) + 2] [2 In 2 - In \Or{t)\] ~ fi^m In \P{t)\ - [(3{t) + 2] In \(3{t) + 2| « 0. 
Notice f3{t) depends on fij, and Or but not on nt- For Or ^ 0{e) ^ 1, /3(i) can be approximated by 



Mb + 1 



(21) 



(22) 



which is between and -1 for /if, > 1. For details of the derivation of (191 and (20), as well as some 
numerical results regarding (22), see Appendix [A| From (20), one has Uc ~ 0{e^) and Vc, flc ~ '^(e'^) 
since Or, Or ^ 0(e), which is consistent with the observation mentioned before: Uc, Vc, i^c 0(e). Since 
= Oc ^ 0{e^) ^ 1, one has Oc ^ 0(1) for a period of time t ~ 0(1). Hence, the transformation 
matrix Rg^ given in ([!]) is approximately an identity matrix, which means the inertial frame velocity can 
be approximated by the body frame velocity, i.e. gc ~ ^c- One can easily obtain the velocities for other 
combinations of the signs of Or and Or by symmetry. 

As an example, when Or is given by (18) and 0„iax = e = 0.1, the velocity of C is numerically solved 
from the full nonlinear equations (11), and the result is plotted in Figure |6] The large- amplitude solutions 
plots (Figurejsj) approximately contain those in Figure|6] Specifically, Or for Omax — 0.1 nearly coincides with 
the two subintervals where Or is between [—0.1 0.1] for the case with 0„iax = 7''/2, when time is dilated by a 
constant factor. Hence, the results in Figure [6] also nearly coincide with the portions around zero in Figure[5] 
One can observe that the order of magnitude of gc is much smaller than e = 0.1. For the first quarter period, 
ic oc (1 — it), T/c cx (1 — 4t)^ and 0c oc (1 — 4t)^, which is consistent with the analytical solutions in ( [20| . 
When t is close to 1/4 and 3/4, the velocities can no longer be approximated by linear or quadratic functions. 
This is because /3{t) can no longer be approximated by a constant since Or is discontinuous, but the velocities 
are still bounded. Due to the symmetry in Or in the four quarter periods, the distance traveled by C and 
the total work during one period are given by 



Xr dt ; 



^J■b 



1 



W 



-f • ^/i„dsdt w ^^te. 



since i/c dt ; 



8int + 1) 

0. Therefore, the efficiency for small-amplitude actuation is approximately 



(23) 



(24) 
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0.5 1 0.5 1 0.5 1 

(a) Xc vs. t (b) j/c vs. t (c) 9c vs. t 



Figure 6: Inertial frame velocities of the center of mass C for the two-Unk model with 6r given by |18| and small amplitude 
of actuation Omax = 0.1. The coefficients of friction are /ij, = 1.3 and nt = 1.7. 



This shows that the efficiency is maximized when fit and /.tb — > oo (i.e. {fib — + 1) ^ !)■ Note 

that the linearization is only valid when fit ^ O(e^). When fit is comparable to e^, which means the body 
is almost frictionless in the transverse direction, the terms in W which are of higher order in e cannot be 
ignored, and these higher-order terms depend on fih as well. This means that transverse friction provides 
the leading-order contribution to the work, and tangential (forward and backward) friction is comparable at 
higher order (see Appendix [A| for details). However, (24) is valid for e — ^ with fixed fit, which indicates 
e — >■ 0, or small amplitude actuation is energetically inefficient, and the efficiency increases linearly with e. 
Thus, maximum efficiency occurs at large amplitude, where geometric nonlinearities play a role. 



Large amplitude optimization We now investigate the case when the amplitude of actuation is not small 
in general. Consider a periodic actuation of the relative angle that varies in the interval 9r G [Omin,Qmax\ 
during t G [0, 1]. Without loss of generality, assume ^r(O) = 0r(l) = Omax, and 9r{tmin) = 0min- To simplify 
the analysis, we consider the case that O^ax and 6min are the only extrema of 9r during the period. In 
other words, 0r varies monotonically between the maximum and minimum, i.e. Or < when < t < imm, 
and 9r > when imm < t < 1. In general, if more local extrema exist, one can always divide the period 
into multiple parts at the extrema, and the following analysis can be modified accordingly. Recall that our 



system has the modified kinematic reconstruction equation (13). That is, for the two-link problem. 




V*{0r,Sr) I 0. = A( 



. Sr 



(25) 



where the exact forms of the components [/*, V* and f2* can be derived from the equations of motion (11 ) 



We will show mathematically that the trajectory of C only depends on the path of 9^ but not on the speed 
9r along the path. For a prescribed 9r, one has 



rt pt pt i-e^it) 

,(t) = / 0^{i)di^ / nc(t)di= / n*{9r,Sr)9rdi^ / 

Jo Jo Jo 

r{t). 



, Sr ) d9r 



-l)d0„ 



for <t <tr 



n* {9r, -1) d9r + Jg"-^'^ n* {9r, 1) d9r, for i„,„ < t < 1, 



(26) 



where t and 9r are integration variables. Once 9c{t) is obtained, the position of C can be given by (14) 



Xc{t) 

ydt) 



ReAcdt- 



0At) / 



cos Uc 

\ sin 9c 



f cos 9 
Vsin6', 

— sin 9, 
cos 9c 



— smf 
cos 9, 



U* (9r, Sr) 
V*{9r,Sr) 



dt 



d9r 



COS Vc 

sin 9r 



— sin 9^ 

COS 9r 



U*{9r,Sr) 
V*{9r,Sr) 



9r dt 



(27) 
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where the integration can be evaluated similarly to (26). Therefore, gc{t) only depends on the path of 9r 
via the integration limits, and the speed 9r does not explicitly appear in the expression (although its sign 
does appear). For distance d, work W and consequently efficiency e, the integration is over the whole period 
from i = to 1, and hence they only depend on the extrema of Or during the period: 9max and 0mi7i- This is 
because the parameter space of shapes is only one dimensional, so the path of Or is defined by the endpoints 

Omax and Oraiji- 



Due to the nonlinearity of the system, the calculation of gc(i) via ( 26 ) and ( 27 ) is not trivial. However, 
from the numerical results of Omax = 0.1 and 7r/2, one can observe that the orientation of the snake Oc 
is usually very small during the locomotion. This is due to the symmetry in shape about the by axis, 
and the fact that the /if,-/i/ asymmetry has little effect on rotation. Indeed, even for a large amplitude 
Omax — TT — 0.01, ouc Still has supt ll^dl < 1-4 degrees for all time. Hence, one can closely approximate the 
problem by assuming ~ during the period. Therefore, from Q, the velocities of C in both frames are 
approximately the same, i.e. gc ~ ^c, in which the non-zero components are the linear velocity, 



U* {Or, Sr) 
V*{Or,Sr) 



(28) 



Since the Heaviside function appears in the tangential but not the transverse direction, comparing the force 



equation in (111 for the same Or and opposite Or gives 



U* 



1) = ~U*iOr, -1), V*{Or, 1) = V*iOr, -1). 



(29) 



This is also evident from the results shown in Figure [5]( a) and (b): for two instants symmetric about 0.5, 
which have equal Or but opposite Or, the corresponding Uc are nearly equal but Vc are nearly opposite. To 
simplify the notation, denote u{Or) = U*{Or,Sr — 1), and v{Or) = V*{Or,Sr — 1). They are nonlinear 
functions of Or only. From (27) and ( [28| , integrating the velocity in the inertial frame for the whole period 
yields 



Xcil) 




U*{0r,Sr)Or At ^ 



-U{0r)0r dt 



U{0r)0r dt 



u{Or) dOr — 2X{9max) — 2X{0min), 



^ V*{0r,Sr)0rdt^ [ v{0r)0rdt+f v{Or) Or dt 



(30) 



V{0r) dOr = 0, 

where X(0) = u{0) dO depends on the form of u and the integration limit 0. The distance d is given by 

d = \/xe(l)2+y,(l)2 « 2X{0max) ~ 2X(0„„„), (31) 

since Xc{0) — yc(0) = 0. The power, or rate of work, done by the snake at a given time can be written 



P{9r,0r)= / -f-iunds^p*{9r,Sr 



(32) 



Power is identical for the instants with the same value of Or and same magnitude but opposite sign of Or- 
That is, P{0r,0r) — P{0r,—0r)- Thcrcforc, one has p*[OrA) = —p*[Or,—l)- Similarly to our definition of 
u{9r), we define p{Or) = p*{Or, 1). The total work is given by integrating power over the period. 



W 



P{0r,0r)dt 



p{Or) dOr = 2W{0„ 



(33) 
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(a) X vs. 9max 
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Figure 7: Two-link model (a) X, (b) W and (c) e as functions of maximum amplitude of actuation Bmax for fii, = 1.3 and 
^lt = 1.7. 



where W{9) = p(6') d0. The integrals X and W only explicitly depend on their upper integration limits (and 
implicitly on the parameters /if, and fj,t). By definition, they are odd functions, for instance, 1L{9) — — X(— 0). 
For symmetric cases, Omax — —Smin (e.g. the one shown in Figure |3]), and the efficiency is simplified to 



max 

max 



(34) 



Figure [7] shows X,W and e as functions of the maximum amplitude 9max for parameters = 1-3 and 
Ht = 1-7. One can see that X is maximized around Omax ~ 37r/4, but since increases faster at larger ^maa;, 
the efficiency e is maximized around 9max ~ 7r/2. We emphasize that, since the parameter space of shapes 
is only one dimensional, as long as 9max is the same, the efficiency and trajectory traveled by C are the 
same. For example, any function 9r that has the same maximum and minimum (and varies monotonically 
in between) as the one shown in Figure [3] for example a cosine function with period 1 and amplitude 7r/2, 
will also result in the trajectory shown in Figure [5j a) and have the same e. 

We repeat the process above to calculate e as a function of 9 max for various /i;, and ^t, and the 
results are shown in Figure [s] for (a) pn, — 1.3 and various and (b) Ht = 1-7 and various ^b- In (a), the 
efficiency-maximizing 9max ~ 7r/2 for /ij > 1, and increases as jit drops below 1, in which case there is a 
smaller increase in work for the additional transverse motion associated with a larger amplitude. In (b), one 
can see that for /it = 1.7 > 1, the efficiency-maximizing 9max ~ 7r/2 regardless of the value of /it. Note that 
when 9max is small, e varies almost linearly with 9max , which agrees with the small-amplitude result in ( 24 1 . 

We now determine the values of /if, and that maximize e. Figure [9|^c) shows a contour plot of e 
as a function of /if, and /if when 9max = 7r/2. Here e is evaluated at points on a 10 x 10 logarithmic grid 
with nodes spaced by factors of 1.5, and /if, ranges from 1.5" = 1 to 1.5^ ~ 38.44 and fit from 1.5~^ w 0.088 
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Figure 8: Two-link efficiency e as a function of actuation amplitude Bmax for various /t;, and /Jt: (a) fixed /j;, = 1.3 and various 
/xt and (b) fixed fit = 1.7 and various /ij. 
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Figure 9: Contour maps of (a) distance d, (b) work W and (c) efficiency e as functions of /ii, and /it (on a log-log scale), for 



to 1.5^ w 3.375. The highest efficiency occurs at the largest /z;, and the smallest /zt in this range. In 
Figure [9j a), the distance d increases when and increase. By contrast, work W is mostly independent 
of and increases when nt increases, as shown in (b). Note this is consistent with the small amplitude 
result in (23). The trends shown in Figure |9] hold for a wider range of parameters, which suggests that the 
efficiency-maximizing parameters are jiy, large and jit small. 



4 Three-link model 




Figure 10: Three-link snake model; see text for descripti( 



We now consider a snake model with three links (each with length 1/3) connected by hinge joints as 
depicted in Figure [TOj The shape is now given by two relative angles, 9^1 and 9^2- Since the links cannot 
penetrate each other, the angles are constrained to lie in the set 



ri,6'r2) e (-7r,7r) X (-7r,7r)P| 



/ 61^2 > -TT - 61^1/2, -TT < 9rl < -27r/3 \ ' 

9r2 > -2tt -29ri, -27r/3 < 0^1 < ~7r/2 

0^2 < 27r + 26*^1, 7r/2 < 6*^1 < -27r/3 

\ 9r2 <TT- 9rl/2, 2tt/'S < 9rl < TT J 



(35) 



In Figure 11 the infeasible regions are shaded at the upper right and lower left corners in the {9ri , 9r2) 
plane. One can describe the motion in the incrtial frame {e^jGy} or in the body-fixed frame {bjjjbj^}, with 
the angle between b^; and now given by 



1 



1 



1 



7r2- 



(36) 



Here 9t, 9m and 9h are the orientations of the tail, middle and head links in the inertial frame, respectively. 



11 




Fi gurG 111 Shape change parameter plane {0r\^0r2) ■ The shaded areas are infeasible due to the mutual avoidance of the links. 
A periodic kinematics is a directional closed path r), with the initial state marked by o. 



In the body frame, the configuration variable on each of the three hnks is 



\s ~ 5/18) cos Gt - l/6cose,„ ~ 1/18 cos 6^^ 
Gf = I (s -5/18) sin Gt- 1/6 sine™ -1/18 sin Bft 

^1/18 cos Gt + (s - 1/2) cos e„ - 1/18 cos e?,^ 
1/18 sin Gt + (s- l/2)sine,„ - 1/18 sin 9,, 

6„, = 0rl/3 — f?r2/3 



(37) 



^l/18cosef + l/6cose™ + (s - 13/18) cosO/i^ 
1/18 sin Gt + 1/6 sin 9™ + (s - 13/18) sin 9^ 

9/i = 0rl/3 + 20J.2/3 

The hnear and angular velocities in the body and inertial frames and the transformations between them 
are again given by (|4| and ([5|. The expressions for forces, the equations of motion and the definitions of 
distance, work and efficiency are all in the same as in the two-link case. 

In this section, to keep the parameter space tractable, we set /i;, = 1.3 and [it — 1.7 (based on the 
experimental measurement in |13| ) and search for the efficiency-maximizing shape change in terms of Q^i (t) 
and 6r2{t)- Now the path of the shape change over one period in {6ri , 0r2) space is a directional closed curve 



(path) denoted by 77 (an example of which is shown in Figure 11) 



The reconstruction equation ( 13 ) is now 



Ui{9rl, Srl,0r2, Sr2) C^2 (^rl, Sri, Qr2, Sr2) 
V]*(6'rl, 5*1.1, ^r2j 3^2) V'2*(6'rl, Sri, (^r2, Sr2) 
{9rl, Srl,0r2, Sr2) ^li^rl, Srl,0r2, Sr2) 




A.(^0r, Sr'jOr- 



(38) 



As in the two-link case, most quantities of interest are independent of how time is parametrized. For example. 



.it) 



1 \"rl, Sri, (^r2, Sr2)(^rl + ^^.^firl, Srl,0r2, 'S'r2)^r2 
/ [^^1 {drl, Sri, (^r2, Sr2) dOrl + ^2i^rl, Srl,dr2, Sr2) d0r2] , 



At 



(39) 



where rf^ is the portion of the shape-change path connecting (0ri(O), 0r2(O)) to {Ori{t),9r2{t)). Hence, 9c{t) 
does not depend on the speeds of the shape change, Ori and 9r2- Similar results hold for Xc{t) and ydt), and 
consequently also for the distance d, work W and efficiency e. Hence, e is only a function of the path 77. In 
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contrast to the two-link case, the shape of the three-hnk model is not symmetric in general. Therefore, the 



orientation 9c given by ( 39 1 is no longer always small, and the initial and final orientations in one period are 
not necessarily the same. 

The extrema of e occur when the variation of e with respect to ?/ is zero, i.e. 



5rie = 0. 



(40) 



Our goal is to find the rj which are local and/or global maximizers of e. 



Triangular waves for Ori and Or2 We start by maximizing within low-dimensional subspaces of the space 
of all feasible paths. First, we prescribe Ori and dr2 to be triangular waves with period 1 (similar to the 
two-link case), and equal amplitudes: Ori^^^ = 0r2mai — ~^rimin — ~^>-2„i„ = (^max- There is a phase delay 
(j) between the two angles: 9ri is assumed to be given by Figure |3][a) and 0^2 is shifted tf> behind Ori in time. 
In general, < 9max < and < < 1, but the constraint given in (35) also applies. The angular velocities 
are given by 



0<t< 0.5, 
0.5 <t< 1, 



-'idmax, O<t-0<O.5or -l<i-0<-O.5, 
46'^„,, 0.5 <t - 6 <1 or -0.5 <t- 6 <0. 



(41) 



Such a path 77 is a rectangle in {9ri, 0r2), centered at the origin with the edges at ±45 degrees to the 9ri axis, 
examples of which are shown in Figure [l2|b) . The shapes and orientations of the rectangular paths lie in a 
two-dimensional space parametrized by {Omax^ 4')- We perform an exhaustive search for the globally optimal 
shape change in this case: discretize {9max, </*) on a 100 x 100 mesh on the range [0.01 , tt — 0.01] x [0.01 , 0.99], 
and calculate the efficiency at each node by solving the equations of motion given in (11). The results are 
depicted as a contour plot of e as a function of Omax and 4> in Figure 12 'a). The areas bounded by the dashed 
lines and the boundaries at the upper and lower right corners are infeasible due to the constraint given in ( 35 ) . 
One finds 3 local maxima of e, located at pi,P2 and p^. At pi, 9,nax — 2.4694, = 0.1981, and e — 0.1197, 
which is the global maximum for this family of kinematics. The second and third local minima are Cp^ = 



0.0834 at 9^ax = 1.6181,0 = 0.5940 and e^^ 
trajectories 77 in the {9ri , 9r2) plane are shown in Figure 
other two are clockwise, because 4> < 0.5 for pi and (j> > 0.5 for the others. Figure [13| shows trajectories of 



0.0634 at 9raax = 2.7847,0 = 0.8118. The corresponding 
h). Note that r/p-^ is counterclockwise, and the 
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C and snapshots of the snake at these local optima. The arrows indicate the directions of travel along the 
trajectories. We note that on average the snake moves forward for pi and p2, but backward for p^ (keep in 
mind that distance is the norm of net displacement). It seems suboptimal to move backwards (since /ifc > 1), 
and, indeed, Cp, is not a global maximum. One can see two small "loops" in the trajectory of C for pi and 




(b) paths in {dr\,dr2) space 



Fi gure 12: (a) Contour plot of e as function of 9max and (f) for the shape change prescribed in ( |41[ l . Areas bounded by dashed 
lines and boundaries at the upper and lower right corners are infeasible due to the mutual avoidance of the links. The three local 



maxima of e are marked pi,P2 and ps: Sp-^ 
paths 77 in the (6rl,^r2) parameter space. 



0.1197 (global maximum) 



0.0834 and e„ 



0.0634. (b) Corresponding 



13 




Figure 13: Trajectories of C (left) and snapshots (right) in the inertia! frame for pi,P2 and p^. Starting at (0,0) at t = 0, 
the positions of C at t = 1 are: for pi, (0.1649, 0.0504); for p2, (0.0536, 0.0105); for pa, (-0.0973, -0.0300). The snapshots are 
given at time increments of 0.05 for pi and p3 and 0.1 for p2. The head of the snake is represented by O. The orientations of 
the snapshots in the inertia! frame are preserved, whiie the centers are shifted to a straight line for better illustration. 



they correspond to the first and third edges of the rectangle, during which the head and tail links are almost 
parallel to each other during the motion. The other two edges correspond to the head and tail links rotating 
in opposite directions, which results in much greater displacements compared to the small loops. During 
the period, the orientation of the middle link does not vary as much as the other two links. The differences 
among the edges in the work done are not very large but the differences in distance traveled are. The first 
and third edges are somewhat analogous to recovery strokes, with the second and fourth edges analogous to 
power strokes. 



More general shape changes Now, we optimize in a more general class of shape changes. Since dri and 
9r2 are periodic functions with period 1, they can be represented as Fourier series, 



n 



^it)^f+Y, [a} COS i2jnt) + h] sin {2jTit)\ , (42) 



where a* and 6* are Fourier coefficients, their upper-scripts i ~ 1 ov i ~ 2 correspond to 9r-i or 0,-2, 
respectively, and the lower-scripts j indicate the mode in Fourier series. One needs n — ^ cxj to represent a 
general function. Following |25| , we start with n — 1 and increase n systematically, observing the resulting 
change in the optima. 

For n — 1, the relative angles are given by Ori{t) = a\/2 + a\ cos (27ri) -I- h\ snv{2nt) and 0r2{t) = 
al/2 + a\ cos {2nt) -f- b\ sin {2nt). The corresponding 77 is a relatively smooth curve, an example of which is 



shown in Figure 14 ^a). For compactness, we denote q = [a\ , , a\ , a\ , b\ , as the coefficient vector. 



Mathematically, the optimization problem can be stated as the following, 

Objective max e(ry), 

q 

Variables q = [aj , Oq , aj , , 6} , , 

77 : Oriit) = ay 2 -f a\ cos (27rt) + b\ sin (27r<) , i = 1, 2, (43) 

Equations /q^ f ds = 0, X-^ • f ds = 0, Vt G [0 , 1], 

Constraints (dri,0r2) & Sg 
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(a) path ri in {9^1, 9r2) space (b) trajectory of C 

I II III IV V 

(c) snapshots 

Fi gure 14: (a) Global optimal path 77 in {9r\, ^r2) space for 1 mode in the Fourier series, (b) Trajectory of the center of mass 
C in the inertial frame, which moves from (0, 0) at t = to (0.3381, —0.0780) at t = 1, with an arrow showing the direction of 
locomotion, (c) Snapshots of three-link snake are given at time increments of 0.05, with the head of snake represented by O. 



where Sg^^.e^^ is given by (35 1. Since q has 6 degrees of freedom, performing an exhaustive search is 
more expensive now. Instead, we utilize the subroutines provided by the Global Optimization Toolbox 
in MATLAB. The procedure is briefly described here. For a random initial guess qo, a local constrained 
optimization function fmincon, which implements a Sequential Quadratic Programming (SQP) algorithm, 
is called to solve for a local optimal q;oca; that results in a local maximum of efficiency eiocai- To find 
the globally optimal variable qgiobai, two approaches are taken: (i) repeat the local search 2000 times with 
random initial guesses pre-filtered such that the new searches do not fall back to the immediate vicinity of 
exploited results; then pick the largest eiocai as a candidate for the global maximum (using MultiStart); 
(ii) numerically calculate the basins of attraction of the local maxima and then pick the largest eiocai as 
a candidate when the basins cover the variable space (using GlobalSearch). Each approach is repeated 
several times, and the largest candidate is regarded as the global maximum. We verify the global maximum 
is indeed a local maximum by numerically computing the gradient using PatternSearch. 

The globally optimal coefficients are qgiobai = [-0.1635, -0.2274, -2.0246, -0.1079, 2.2364, -2.9886]^, 
and the corresponding 77 is depicted in Figure [T4{ ^a), with the starting point marked by o and the direction 
shown by the arrow. Figure |14[ b) shows the trajectory of C in the inertial frame, which moves from the 
origin when i = to (0.3381, —0.0780) when t ~ 1. Figure 14 |c) shows the snapshots of the snake at time 
increment At = 0.05 from i = to 1. The five instants t — 0, 0.25, 0.5, 0.75, 1 correspond to the five states / 
to V, respectively. The distance is c? = 0.3470, the work is 14^ = 1.5123, and the efficiency is e = 0.2295. In 
comparison to the optimal solution in the triangular wave case (e = 0.1197), the efficiency is much higher. 
The loops in the trajectory in Figure [T3| no longer exist, the distance is much larger, and the trajectory is 
smoother than in the triangular wave case. 

For n = 2, the coefficient variable is q = [aj , aQ , a\ , of , bl , bf , , , b2 , b'^]'^ , and the 4 additional 
coefficients correspond to the second mode in the Fourier series. Following the aforementioned procedure, we 
obtain the global optimal coefficient q = [-2.7831, -2.0203, 0.3941, 0.0653, -0.3580, 0.4929, -0.0156, 0.1473, 
0.1589,-0.0305]"^. Note that the second-mode coefficients are smaller on average than the ffi-st-mode coef- 
ficients. Besides this global maximum, we also find 3 other local maxima of interest. Their paths in the 
parameter space are plotted in Figure 15 a). We denote the global optimal path at the lower left corner 
as 771, and the local optimal results at the upper right, lower right and upper left corners as 772, 773 and 774, 
respectively. The corresponding efficiencies are ei = 0.3253, 62 = 0.3252, 63 = 0.2893 and 64 = 0.2882. Note 
that they are all larger than the global maximum with one mode. One can see that 771 and 7/2 are very similar 
in shape and are counterclockwise. Their efficiencies are almost equal. Similar results hold for r/3 and 774, 
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(a) paths in parameter space 



(b) trajectories and snapshots of rji and 773 



Fi gure 15: (a) Optimal paths rji^i = 1,2,3,4 in (^?rii^r2) space for n = 2 mode in the Fourier series; ei = 0.3253 (global 
optimum), 62 = 0.3252, 63 = 0.2893 and 64 = 0.2882. (b) Trajectories of the centers of mass C in the inertial frames for jji (top), 
for which C moves from (0, 0) at t = to (0.0884, 0.0144) at t = 1. The snapshots are at t = 0, 0.185, 0.37, 0.555, 0.74, 0.87 and 1. 
For ri3 (bottom), C moves from (0, 0) at t = to (0.0586, 0.0454) at i = 1. The snapshots are at t = 0, 0.19, 0.38, 0.555, 0.73, 0.875 
and 1. 

only they are clockwise. In all four cases, the path shapes are close to 45-45-90 right triangles. The locations 
of the centers of 771 and r/2 are close to the 9ri = 9r2 hne and those of 773 and 774 are close to the 6ri = —9r2 
line. Figure [l5|b) shows the trajectories of C in the inertial frame and snapshots of rji and 773. For 772 and 
7/4, they are nearly mirror images of those for 771 and 773, respectively. In all four cases the motion is the 
following: moving just the head or just the tail link first, then moving both, then just the other link. The 
orientation of the middle link does not change much in all cases. It is interesting to note that for the 2-mode 
case, more efficient kinematics correspond to a "contraction-expansion" motion, which is reminiscent of the 
"concertina" mode of snake locomotion |17j . 

For 71 > 2, so far we have not found any motion that is more efficient that the ones found for n — 2. 
The numerical error in our calculation of efficiency is less than 10~^, based on the time discretization used. 
The higher dimensionality of the parameter space at larger n implies more computational time is needed to 
locate local optima, which may also increase in number. Nevertheless, based on our computations for tt, > 2, 
the global optimum obtained in the 2-mode study seems likely to be close to the global efficiency-maximizing 
shape change for the three-link model. 

5 Discussion 

To summarize, we have adapted a model for the slithering locomotion of snakes from [T^] to the case of 
two-link and three-link bodies sliding in 2D. Two dimensionless numbers, the ratios of the coefficients of 
friction, are the key physical parameters. Because of the frictional anisotropy, the local connection matrix 
in the reconstruction equation depends on both the relative angles and the directions of their velocities, 
which breaks the symmetry of time reversal and shape reversal equivalence, and the Scallop theorem does 
not apply. We maximized the efficiency e, the ratio of the distance traveled to the work done by the snake 
during one period of actuation, under various kinematic assumptions. For a two-link snake in the limit 
of small-amplitude actuations, e is maximized when /x;, ^ 1 and <C 1. Simulations of large-amplitude 
actuations show that the optimal shape change occurs when the relative angle amplitude Omax ~ 7r/2 except 
for small /Xj. In real snakes, /ij is usually not very small. We note that when /Xj <C 1, the distance traveled is 
also very small, as shown in Figure [9] Although the body can increase its actuation frequency by decreasing 
T to achieve higher d in a given time, the actuation frequency may be constrained in biological or robotic 
snakes. If so, a very small may constrain the speed of locomotion to be small, which may be undesirable 
for biological or robotic snakes. We also note that Hb can be increased by real snakes by altering the angles at 
which their scales contact the ground |17j . For the three-link model, we searched for optimal motions in terms 
of the two relative angles O^i and 6^2, a-nd assumed = 1.3, /ij — 1.7, which is a set of coefficients measured 
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for juvenile milk snakes [T3]. We first studied the family of kinematics with relative angles constrained to 
a rectangular path in parameter space. Three local optima were found, and the global optimum results in 
a trajectory with sharp changes in center-of-mass velocity, shown in Figure |13[ We then considered paths 
parametrized by Fourier series for 0^1 and 0^2- If only the first mode (lowest frequency) is allowed, the 
optimal shape-change path in the parameter space resembles an ellipse elongated in the diagonal direction, 
shown in Figure |14[ If more modes are allowed, then the optimal paths are instead close to right triangles 
that are small compared to the ellipse, shown in Figure 15 The globally optimal motion is found to be 
reminiscent of the concertina mode of snake locomotion |17j . Table [T] highlights some of our findings. 



Model 



Optimal 
shape change 



Friction 
coefficients used 



Optimal 
friction coefficients 



Two-link 




unless /J,t ^ 1 



1 < Aifc/M/ < oo 
< ^it|^J'f < oo 



Three- link 



8rl 



^ib/^J'f = 1-3 

fit/nf ^ 1.7 



Table 1: Highlights of resuhs for two-link and three-link snake models, dr or 9ri and 9r2 are the angles between the links, 
and fJ.f,fJ.b and fit are the forward, backward and transverse coefficients of friction, respectively. In the two- link model, the 
optimal coefficients are obtained through analysis; in the three-link model, the coefficients of friction are fixed to be those from 
an experiment |13| . 

Although the results obtained in this work are based on a simplified model of snake locomotion, the 
system is still mathematically challenging due to its nonlinear nature. A key feature of snake locomotion 
- frictional anisotropy - is present in this work, but to keep the analysis tractable, other biological and 
dynamical aspects of snake slithering are omitted. One aspect is the dynamic load distribution, i.e. snakes 
lifting part of their body during locomotion [T^]. To take this into consideration, the current model can be 
modified such that coefficients of friction can change in time and along the snake. This modification will 
result in a more complex optimization problem. The assumption of Coulomb friction could also be replaced 
by more complex and nonlinear friction models such as those discussed in |26j , which would also present new 
challenges in solving the dynamical equations as well as the optimization problems. 

In the three-link case, we focused on the optimal shape change with fixed coefficients of friction; a 
natural extension is to vary these. Another natural extension is to systematically increase the number of 
links N and see how the optima change. For large N it would interesting to compare optimal motions with 
those found in a two-parameter space of smooth shapes by [13] . 
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A Derivation of the equations of motion and small-amplitude anal- 
ysis 

The reader is first reminded that the velocity of an arbitrary point on any link can be decomposed into the 
velocity of the link center and a rotation about that center (see Figure [2]). For the tail link, in the body 
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frame, the link center velocities in the tangential and transverse directions are 



JJ7 



''cos 


y 


~sm-\ 




ysin 


y 


COS - ; 





1^ dr 
-lie COS — 

4 2 



(44) 



and the rotation rate is 8f = — 0,./2. The unit tangent vector of the tail link is St = (cos(0r/2) , — sin((?r/2)), 
and the unit transverse vector is fit = (sin(6',./2) , 008(6*^/2)). Therefore the linear velocity of any point on 
the tail link is given by 



Similarly, for the head link, the link center velocity is 

/ 



COS — sm — 
2 2 



\ — sin — COS — 
^22 



fit, for < s < ^. 



1 ■ I 

Ur — -Or sin - 



(45) 



(46) 



-51c cos — 
4 2 



and Qh = ^r /2, the unit vectors Sh — (cos(0r/2) , sin(^?j./2)) and fi/i = (— sin(6',./2) , cos(0r/2)). Therefore 
the velocity of any point on the head link is given by 



The unit velocity vector for the tail link is 



n/i, for ^ < s < 1. 



1 



ur -{s-i) 



4/ 4 



(47) 



(48) 



and that for the head link is of a similar form. For concreteness, we discuss the case when Or > and 
Or < 0, as depicted in Figure [2} The analysis can be easily extended to other cases with minor modifications. 
Substituting ( 45 ) into ^ and integrating over the whole link, the frictional force exerted on the tail link is 



i: 



ids 



Of. 



(ur + Qt/^) + (c/f)2 - ^ ([/,"- 6^/4) +{uir 



In 



up + Qt/^ + Jlui^ + et/A] +{U!r 



(49) 



Similarly, the frictional force on the head link is 



ids 



Ok 
Ok 



uj^ + eh/^ + xliK + e^/A] +{Ufir 



(50) 



The force equation (11) is ids + J^^^tds — 0. One can readily derive the torque equation, which is 
omitted here due to complexity. We have given the force and torque equations in terms of , C/" , and 
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UJ^. The center of mass velocity — [t/c, f^c]"^ can be easily obtained by inverting (44 1 and (46 1. In 
general, does not have a closed- form solution for large amplitudes of 9r- 

For small amplitude kinematics supj = e 1, one can make the approximation 



cos(6i,./2) « 1, sin(6'^/2) w 61^/2. 



(51) 



Hence, Sj and Sh are almost parallel to b^^ while fif and hh are almost parallel to bj,. One can intuitively 
see that J7c, Vc, ^ 0{e) due to the symmetry in kinematics for the two-link model, and this can also be 
verified later. The dominant term in velocity is due to the rotation about each link center. As an example, 
in the tail link (45), the dominant term is (s— l/4)0(nt, whose magnitude is of the order 0(e) and direction 
is almost parallel to hy. All other terms have order of magnitude O(e^), and they are comparable to the 
dominant term only when s ~ 1/4. Therefore, the components of the unit velocity vector in ( |48p are: almost 
sign function in hy , and mostly zero in b^; except around s « 1/4. A similar result holds for the head link. 

We now analyze the orders of magnitudes in the equations of motion. For the b^^ component of the 
force equation, the term generated from the transverse direction, /it6'^/2 ^ O(e^), is negligible compared to 
the tangential direction terms ^ 0(e), as long as fJ-t/fJ-f and fJ^b/l^f are 0(1) in the small-e limit. Hence, after 
multiplying by the common denominator 9r/2 and neglecting higher order terms, the equation is reduced to 



Mb ( t^c + Y^drOr 1 In 



Uc - 1 In 



4Ur - 



./4 



(52) 



For the by component of the force equation, although the absolute values of the transverse terms are 0(1), 
they almost cancel out due to the symmetry in shape. As a result, contributions from both tangential and 
transverse directions are comparable in order of magnitude and ~ O(e^). Neglecting higher order terms, the 
equation becomes 



Ait ^V, 



^Mlj - /ib (uA- + ] In 



AUc + 0r0r/4. 



0, 



(53) 



in which (52 1 was used. For the torque equation, after multiplying by the common denominator 0j./2, 
contributions from the tail and head links are comparable and ^ 0{e^). After manipulation, the linearized 
torque equation is reduced to 

3 



-Ait {2Ucer - n^) 



4Mt(l + Aifc)^ + f (l-Aih)^^ 



Uc + —erOr I In 



4Uc + 9r9r/4: 



0, 



(54) 



in which (52 ) was used again. To solve the 3 equations (52 H(54), notice in (52 ) that Uc appears with 9^0^/16, 
and hence we seek Uc in the form 

1 



Uc 



-^[l + /3(t)] 



(55) 



where /3(t) is a function of time. The governing equation for /3(t) can be obtained by substituting (55) 
into (52), obtaining 



Mb 



16 



In 



1 

16 



[2 + I3{t)] Mr In 



[2 + m] 



0, 



(56) 



which can be further simplified to 

[(^ifc + l)l3{t) + 2] [21n2 - In |0,(t)|] - fi^m In \l3{t)\ " + 2] In |/3(t) -t- 2| « 0. (57) 



Note that /3 depends on but not on For 9r ^ 0(e) <C 1, the dominant term in (57) is In 16*^(^)1. Hence, 
for the left side of (57) to be bounded in the small-e limit, the coefiicient of the dominant term has to tend 
to zero, i.e. 

t + Mb 
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m 

-0.5 r 



t 

1/8 1/4 

Fi gUre 16! /3 (solid) as a function of t for ^5 = 1.3 and e = 0.1. It is close to a constant — 2/(/i^ + 1) = —0.8696 (dashed). 



For the parameter /it, — 1.3, and 9r{t) given by the triangular wave in (18) for e =^0.1, when > 

One can see that 



16 



and 6r < (first quarter period), /3(t) is calculated from (57) and plotted in Figure 
— 2/(/ib + 1) = —0.8696 is a good approximation of P{t) for most times. One can obtain solutions of Vc and 
r^c by substituting ( 55 ) into the remaining two equations of motion: 



1 

32 



-1 



(3 In 



-^(l + /3)-^/3[/3(l + A^b)+2/ib] In 



(59) 



As stated before, these solutions for Uc, Vc and flc are for the case when 9^ > and 9^ < 0. One can 
easily obtain solutions for other cases based on symmetry. The reader is reminded that one can approximate 
inertial frame velocity with body frame velocity since ^ 1 for all time, i.e. gc ~ ^c- As an example, for 
/If, — 1.3, /it = 1.7 and e = 0.1, the velocities are given in Figure|6] Since the velocity in the ej^-direction yc 
is almost anti-symmetric about t = 1/2, one has ycdt « 0. Hence, the distance is given by 



Xr dt ; 



1/4 



Ucdtxi 

IQ JO JO 

.1/4 ^ 

"^io -I^(l + '5)(-4e)e(l-4t)d< 



1/4 



-—{i+mAdt 
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The work is given by 



W ' 



1 i-i /-i/i / /•1/2 



Jo 











f • Inn ds dt 



1/2, 



1 



(60) 



(61) 



where 



1 



16(/ih + 1)2 



2/it + ^j,b + ^J■l + jj^tyl + A^fc + Mb (m& + 1)^ 2/ig + /it 



2e2 



ln//& 



Therefore, the efficiency is given by 



One can see that when fit ^ 476^ 



_ /ifc - 1 

2(/ib + l)(/it + 47e2)' 

e ~ re. 

2/*t(A*f) + 1) 



(62) 



(63) 



(64) 
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